#!/usr/bin/perl -w
use strict;

# To use this script you need to have run HMMSEARCH or 
# HMMSCAN with this option
# --domtblout domains.out
# hmmsearch-2.3.2 HMMFILE PROTEINFILE > PROTEIN-vs-HMMFILE.hmmsearch
#
# hmmsearch-2.3.2 Sugar_tr.hmm Cryptoproteins.pep > Cryptoproteins-vs-Sugar_tr.hmmsearch
#
# bp_hmmer_to_table.pl Cryptoproteins-vs-Sugar_tr.hmmsearch > Cryptoproteins-vs-Sugar_tr.hmmsearch.tab
# the tab file is what you provide to this script!
# to run this you do
# perl hmmer3_extract_domains -i Sugar_tr.crypto.domains.out -db Cryptoproteins.pep 
# when it is done you will have a file called Sugar_tr.domains.out.Sugar_tr.domains.fa

use Getopt::Long;
use Bio::SeqIO;
use Bio::Seq;
use Bio::DB::Fasta;

my ($tablefile,$dir,$database);
my $evalue_cutoff = 1e-5;
my $border = 0;
$dir = '.';
GetOptions('i|input:s' => \$tablefile,
	   'db:s'      => \$database,
	   'd|dir:s'   => \$dir,
	   'e|evalue:s' => \$evalue_cutoff,
	   'b|border:i' => \$border,
	   );

$tablefile ||= shift @ARGV;
mkdir($dir) unless -d $dir;
my $dbh = Bio::DB::Fasta->new($database);
open(my $fh => $tablefile) || die $!;
my %seen;
my %track_seq_domain_count;
my @domains;
while(<$fh>) {
    next if (/^\#/);
    chomp;
    my ($query, $targetstart,$targetend,
	$target,$hmmstart, $hmmend,
	$bitscore,$evalue,
	$description) = split(/\s+/,$_);
    
    next if $evalue > $evalue_cutoff;
    push @domains, [ $query, $hmmstart, $hmmend,
		     $target,$targetstart,$targetend,
		     $evalue, $bitscore,
		     $description ];
    $track_seq_domain_count{$query}->{$target}++;
}

my %recount;
my $n = 1;
# sort by domain, then sequence name, then start position in sequence
for my $d ( sort { $a->[0] cmp $b->[0] || 
		       $a->[3] cmp $b->[3] ||
		       $a->[1] <=> $b->[1] 
		   } @domains ) {
    my ($query, $hmmstart, $hmmend,
	$target,$targetstart,$targetend,
	$evalue, $bitscore,
	$description) = @$d;
    my $n = ++$recount{$query}->{$target};
    my $total = $track_seq_domain_count{$query}->{$target};
#    next if $cevalue> $cevalue_cutoff;
	my $domain = $dbh->seq($target,$targetstart-$border,$targetend+$border);
    
    push @{$seen{$query}}, Bio::Seq->new(-seq => $domain,
					 -id  => sprintf("%s.dom%d_%d",
							 $target,$n,$total),
					 -desc=> sprintf("%d..%d %s:%d..%d %s",
							 $targetstart-$border,
							 $targetend+$border,
							 $query,
							 $hmmstart,$hmmend,
							 $evalue));
}

for my $q ( keys %seen ) {
    my $out = Bio::SeqIO->new(-format => 'fasta',
			      -file   => ">$dir/$tablefile.$q.domains.fa");
    
    $out->write_seq ( @{$seen{$q}} );
}
	
